/*!  \file  BdsDC.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      06/30/2013 02:00:55 PM
 *   \copyright GNU Public License.
 */
#include "BdsDC.h"
#include "Matrix.h"
extern "C" {
void dbdsdc_(
char   *  uplo,
char   * compq,
int    *     n,
double *     d, 
double *     e, 
double *     u, 
int    *   ldu, 
double *    vt,
int    *  ldvt,
double *     q,
int    *    iq, 
double *  work, 
int    * iwork, // *  IWORK   (workspace) INTEGER array, dimension (8*N)
int    * info);

int ilaenv_( const int* ispec, const char* name, const char* opts, 
                 const int* n1, const int* n2, const int* n3,
                 const int* n4 );
}
BdsDC::BdsDC() :BiSVD("Lapack dbdsdc") {
	aUpper    = true;
	aUColMajor = true; 
	aVColMajor = true; 
	setCompileTime(__DATE__, __TIME__);

	int ispec = 9; int dummy = 0;
   smalliz = ilaenv_(&ispec, "dbdsdc", " ", &dummy, &dummy, &dummy, &dummy);
}
void BdsDC::compute(const Matrix & B, Matrix & U, Matrix & Sigma, Matrix & VT){
//	cout << this << endl;
	reset(); success = false;
	assert(U.isColMajor() && VT.isColMajor());
	if ( !B.isBiDiagonal() || !B.isSquare()) return;
	success = true;
	char uplo = B.isUpper() ? 'U' : 'L', compq = 'I';
	int n = B.cols(), ldu = U.lda(), ldvt = VT.lda(), info = 0, lwork = -1;
	double * d = Sigma.getModData(),*e=NULL; allocate(e,n); double * save_e = e;

	memcpy(d, B.getData(), n * sizeof(*d));
	memcpy(e, B.getOffDiagonal(), (n-1)* sizeof(*d));

	int tmp1 = (int)(log(n) - log(smalliz+1));
	int ldq  = n * (11 + 2 * smalliz + 8 * tmp1);
	int ldiq = n * (3 + 3 * smalliz);

//	LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
//	LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))

	double * q  = NULL; allocate(q, ldq);
	int    * iq = NULL; allocate(iq, ldiq);
	int * iwork = NULL; allocate(iwork, 8 * n); int * save_iwork = iwork;

	lwork = compq == 'I' ? 30 * n*n + 40 * n : (compq == 'P' ? 6 * n : 4 * n);
	double * work = NULL; allocate(work, lwork); double * save_work = work;

	// dbdsdc_ sometimes modifies the pointer 'e'
	// as a result when I free the memory, an error occures
	// currently I made a work around by saving the pointer 
	// before calling the subroutine
	// in other case simply causes segmentation fault
	dbdsdc_( &uplo, &compq, &n, d, e, U.getModData(), &ldu, 
				VT.getModData(), &ldvt, q, iq, work, &lwork, &info );
	success &= info == 0;
	e = save_e; save_e = 0;
	work = save_work; save_work = 0;
	iwork = save_iwork; save_iwork = 0;
	delocate(q); delocate(iq);
	delocate(e); delocate(iwork); delocate(work); delocate(iwork);
//	cout << "info = " << info << endl;
//	cout << B - U * Sigma * VT << endl;
//	cout << "aide " << this << endl;
}
//REGISTER(BdsDC);
